%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 模拟退火算法
% 简介：模拟退火算法寻找最大值
% 作者：Zhaojiang
% 日期：2023/10/10
% 企鹅：277746470
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;clear;clc;
%% target fucntion
[x,y]=meshgrid(-4:0.1:4);
z=target(x,y);
mesh(x,y,z);
title('Target Function')
xlabel('x');ylabel('y');zlabel('z');
pause(2);
close all;clear;clc;
%% 设定算法参数
rng(1234);
var_num = 2;% 2个求解变量
lower_bound = -4;% 解的最小值（小数后4位）
upper_bound = 4;% 解的最大值（小数后4位）
bound_length = upper_bound-lower_bound;
step_size = 0.5;% 优化步长
max_iteration = 1000;% 最大降温次数
temp_init = 100;% 初始温度
temp = temp_init;% 当前温度
cooling_rate = 0.99;% 设置冷却速率
%% 算法数据初始化
var_value = lower_bound+bound_length*rand(1,var_num);% 变量值(十进制)
% var_value = [0;2];
func_value = -target(var_value(1),var_value(2));% 函数值(最大值问题转换为最小值)
best_var= var_value;% 最佳变量值
best_func_value = func_value;% 最佳函数值
best_iteration = 0;% 最佳迭代次数
all_best_func_value =zeros(1,max_iteration);% 最佳函数值记录
%% 开始迭代求解
for iteration=1:max_iteration
    % 生成新解
    new_var_value = var_value+step_size*(2*rand(size(var_value))-1);
    new_func_value = -target(new_var_value(1),new_var_value(2));% 最大值问题转换为最小值
    % 新旧解差值
    delta = new_func_value - func_value;
    % 判断是否接受新解
    if delta<0 || (rand >= exp(- delta / temp))
        var_value = new_var_value; % 更新当前解
        func_value = new_func_value; % 更新当前解的目标函数值
    end
    % 是否当前最优解
    if func_value < best_func_value
        best_var = var_value;
        best_func_value = func_value;
        best_iteration = iteration; 
    end
    all_best_func_value(iteration) = best_func_value;
    temp = temp*cooling_rate;
end
%% 处理结果
% 打印结果
max_x = best_var(1);
max_y = best_var(2);
max_ans = target(max_x,max_y);
disp 最优迭代次数
disp(best_iteration)
disp 最大值位置
fprintf("x=%f\ny=%f\nz=%f\n",max_x,max_y,max_ans);
% 可视化结果
figure('position',[200 200 1280 480])
subplot(1,2,1);hold on;
[x,y]=meshgrid(-4:0.1:4);
z=target(x,y);
mesh(x,y,z);
title('Target Function')
xlabel('x');ylabel('y');zlabel('z');
view(-37.5,30)
pause(0.1)
scatter3(max_x,max_y,max_ans,'r*');
subplot(1,2,2);hold on;
x = 1:max_iteration;
y1= all_best_func_value;
try
    assert(length(x)==length(y1))
    plot(x,y1,'-','LineWidth',1.2)
    xlabel('x');ylabel('y');
    title('迭代数据')
catch
    disp("plot failed")
    disp(size(x))
    disp(size(y1))
end